#import matplotlib
import numpy as np
import csv
import random
from datetime import datetime
from time import strftime
from sympy import *
import bisect
#import matplotlib.pyplot as plt

'''
revision of test24
to calculate the false nagetive, false positive

'''

def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_data_list(inputCSV):
    sKey = []
    fn = inputCSV
    ra = csv.DictReader(file(fn), dialect="excel")
    
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        for item in ra.fieldnames:
            temp = float(record[item])
            sKey.append(temp)
    sKey = np.array(sKey)
    sKey.shape=(-1,len(ra.fieldnames))
    return sKey

def nullH():
    c = 6000
    clusterPop = [360275, 1684327, 7627173]
    #n = clusterPop[2]
    N = 29535210
    #k = 1.645 * ((c*n*(N-n)/(N*N))**0.5) + c*n/N
    for n in clusterPop:
        print '----------'
        print 'pop = ', n
        k = Symbol('k')
        k = solve((k-c*n/N)/((c*n*(N-n)/(N*N))**0.5) - 1.645, k)
        k = k[0]
        r = Symbol('r')
        eqn = Eq((((N - n + n * r) * k - c * n * r) ** 2)/(c * n * r * (N - n)), 3.09 ** 2)
        r = 0
        temp = solve(eqn)
        for i in temp:
            if i > 1:
                r = i
        if r > 0:
            print 'r = ', r
        else:
            print 'error in solving r'

        print 'E(c|Ha) = ', c * n * r / (N - n + n * r)
        print 'Var(c|Ha) = ', c * n * r * (N - n) / ((N - n + n * r) ** 2)
        print 'Var(c|Ha) = ', c * n * r * (N - n + n * r - c * r) / ((N - n + n * r) ** 2)

class WeightedRandomGenerator(object):
    def __init__(self, weights):
        self.totals = []
        running_total = 0

        for w in weights:
            running_total += w
            self.totals.append(running_total)

    def next(self):
        rnd = random.random() * self.totals[-1]
        return bisect.bisect_right(self.totals, rnd)

    def __call__(self):
        return self.next()

#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    print "begin at " + getCurTime()

    # total pop: 29535210

    
    # old risk area
    mixed = [91,98,101,104,114,115,119,126,131,142,146,147,154,162,168,172] # pop: 1684327
    rural = [8,9,10,11,12,13,14,15,17,19,20,26,28,33,34,37] # pop: 360275
    urban = [105,107,112,120,122,125,127,128,130,133,134,141,143,149,152,155]   # pop: 7627173
    '''
    # new 2 risk area
    mixed = [95,99,109,110,114,115,117,119,126,131,139,140,142,146,147,237] # pop: 1611198
    rural = [9,10,12,13,14,17,20,23,26,27,28,33,34,35,38,41]    # pop: 501040  
    urban = [125,127,128,130,133,134,136,138,141,143,149,151,152,156,157,163] # pop: 7025156
    '''
    riskarea = mixed + rural + urban
    #print riskarea
    
    #unitCSV = "C:/TP1000_1m.csv"
    unitCSV = 'C:/_DATA/CancerData/SatScan/mult6000/three16_format.csv'
    dataMatrix = build_data_list(unitCSV)  # [id, pop, cancer1, cancer2, cancer3]
    id = dataMatrix[:,0]
    pop = dataMatrix[:,1]
    totalPop = np.sum(pop)
    

    NONriskarea = []
    for i in id:
        if int(i) not in riskarea:
            NONriskarea.append(int(i))
    
    output = []
    output.append(id)
    output.append(pop)
    #print output
    

    x = WeightedRandomGenerator(pop)
    #print x.next()
    # average, mixed, rural, urban

    # old risk setting
    totalCase = 6000
    riskSetting = [0.00020314736208071655, 0.000257364319769158, 0.000334581315540396, 0.000224573495496856]

    #totalCase = 5000
    #riskSetting = [0.000169289468401, 0.000219113454929746, 0.000288423125716719, 0.000188909820040268]

    #totalCase = 4000
    #riskSetting = [0.00013543157472, 0.000180428707141221, 0.000244465947598914, 0.000152920254928375]

    #totalCase = 3000
    #riskSetting = [0.00010157368104, 0.000141050779010198, 0.000199082933979409, 0.000116785036523014]

    # new 2 risk setting
    #totalCase = 6000
    #riskSetting = [0.00020314736208071655, 0.000258669757963189, 0.000310308241128273, 0.000225874034797897]

    #totalCase = 5000
    #riskSetting = [0.00016928946840059713, 0.000219989166072307, 0.000268069531147527, 0.000190056891593152]

    #totalCase = 4000
    #riskSetting = [0.00013543157472047768, 0.000181534121757118, 0.000224937105979011, 0.000154033408421923]

    #totalCase = 3000
    #riskSetting = [0.00010157368104, 0.000141689614757690, 0.000180815027130651, 0.000117713947353513]
    '''
    output = np.zeros(len(id))

    i = 0
    while i < totalPop/10:
        tempID = x.next()
        output[tempID] += 1
        i += 1
    '''
    for i in range(0, 1000):
        print i
        cancer = np.zeros(len(id))
        j = 0
        nonCount = 0
        while j < totalCase:
            tempID = x.next()
            if tempID in NONriskarea:
                if nonCount < 2:
                    risk = riskSetting[0]                
                    nonCount += 1
                    if np.random.random() < risk * 50:
                        cancer[tempID] += 1
                        j += 1
            else:
                if nonCount > 1:
                    if tempID in mixed:
                        risk = riskSetting[1]
                    elif tempID in rural:
                        risk = riskSetting[2]
                    elif tempID in urban:
                        risk = riskSetting[3]
                    else:
                        print 'error a'
                        break
                    if np.random.random() < risk * 50:
                        cancer[tempID] += 1
                        j += 1
                    nonCount = 0
                
            #if np.random.random() < risk * 50:
                    #cancer[tempID] += 1
                    #j += 1
        output.append(cancer)
        #print cancer
    
    #matplotlib.pyplot.hist(cancer, 50, normed=1, facecolor='green', alpha=0.75)
    #matplotlib.pyplot.show()
    output = np.array(output)
    output.shape = (-1, len(id))
    output = np.transpose(output)
    
    #y = mlab.normpdf(50, 100, 15)

    #fileLoc = 'C:/_DATA/CancerData/SatScan/own/'+str(totalCase)+'_modified.csv'
    fileLoc = 'C:/_DATA/CancerData/SatScan/own/temp_6000.csv'
    np.savetxt(fileLoc, output, delimiter=',', fmt = '%10.5f')
    #print output
    #unit_attri = np.zeros((len(dataMatrix),4))
    
    # unit_attri [cancer, pop, rate, id]
    #unit_attri[:,1] = dataMatrix[:,1]
    #unit_attri[:,0] = dataMatrix[:,2]

    '''

    totalPop = 0
    for item in unit_attri[:,1]:
        totalPop += item

    pop = 0
    for item in urban:
        pop += unit_attri[item,1]

    cancerT = []
    for c in range(0,1000):
        cancer = 0
        for item in urban:
            cancer += dataMatrix[item,2 + c]
        cancerT.append(cancer)
    print np.var(cancerT)
    '''
    #print pop*6000*1.15/(totalPop-pop+pop*1.15)

    #nullH()
    
    print "end at " + getCurTime()
    print "========================================================================"  

           
